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Abstract 



Schrodinger equations with nonlinearities concentrated in some regions of 
space are good models of various physical situations and have interesting 
mathematical properties. We show that in the semiclassical limit it is pos- 
sible to separate the relevant degrees of freedom by noticing that in the re- 
gions where the nonlinearities are effective all states are suppressed but the 
metastable ones (resonances). In this way the description of the nonlinear 
regions is reduced to ordinary differential equations weakly coupled to stan- 
dard Schrodinger equations valid in the linear regions. The idea is illustrated 
through the study of a prototype equation recently proposed for resonant 
tunneling of electrons through a double barrier heterostructure and for which 
nonlinear oscillations have been numerically predicted. 
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I. INTRODUCTION 

In a recent paper [2] a nonlinear Schrodinger equation was proposed to model the effect 
of the accumulation of electric charge in the resonant tunneling of electrons through a dou- 
ble barrier heterostructure [3]. The main feature of this equation is that the nonlinearity is 
concentrated only in the region where the resonant state is localised, that is within the two 
barriers. The equation of Ref. [2] is a prototype of a class of nonlinear Schrodinger equations 
which can be used to model several physical situations common in mesoscopic physics, e.g. 
superlattices. In this paper we show that as far as the evolution of the solutions over a 
finite time interval is concerned, it is possible to carry out an adiabatic theory which sim- 
plifies considerably the problem. In fact in the semiclassical limit it is possible to separate 
the relevant degrees of freedom by noticing that in the regions where the nonlinearities are 
effective all states are suppressed but the metastable ones (resonances). In this way the 
problem can be reformulated in terms of ordinary differential equations describing the evo- 
lution in the nonlinear regions weakly coupled to standard Schrodinger equations describing 
the propagation in the linear regions. We shall illustrate this idea, which we think is novel, 
in the prototype case of resonant tunneling through a double barrier. Our approach, beside 
providing an important mathematical advance, allows a better understanding of the physics 
involved. In particular we explain why the oscillations discovered in [2] disappear in the 
limit of very small and very large spacial width of the incoming cloud of electrons while they 
are maximal in intermediate situations. From the mathematical standpoint our discussion 
will be informal and computer assisted. However we think that a substantial part of it can 
be made completely rigorous. It is interesting to remark that the computer time needed 
for the numerical solution of the problem is reduced by several order of magnitudes with 
excellent agreement with the exact results. 

Now we briefly recall the physical origin of our model problem. Interaction among the 
electrons can play a crucial role in electrical transport properties of mesoscopic systems [4]. 
As an example, let us consider a cloud of electrons moving in a double barrier heterostructure 
in which the well region confined between two potential barriers acts like a capacitor whose 
energy changes according to the electron charge trapped in it. We emphasize that the 
localization of the interaction is justified due to the existence of a resonance state which 
permits a long sojourn time inside the double barrier and therefore an accumulation of 
charge. 

The exact time evolution of the depicted system can be written in terms of the anticom- 
muting field operator associated to the electron cloud 

*(r,t) (1) 

where the external potential V(x) = Vo {^[a,b] + l[c,d]J with a < b < c < d depends only on 
the coordinate x orthogonal to the interfaces. As argued in [2] we can assume a decoupling 
between the degrees of freedom parallel and orthogonal to the interfaces and set up an 
approximate description in terms of a Hartree equation for a one particle wave function 
which depends only on the orthogonal coordinate. More precisely we factorize the single 
particle mean field in the following way 
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dt v ' ' 



^V 2 + V(x) + 
Am 



e\r 



,t)V(r',t) df 



V(r,t) = (V(r,t)) ~ u(x,t)(f)(y,z)e-* tj H t . (2) 

where (. . .} means expectation in the many body state at time t; <f)(y, z) is a solution of the 
free particle Schrodinger equation in the plane parallel to the interfaces normalized to unity. 
We consider the interaction term effective only inside the well [6, c] and we approximate it 
by the expression 



— \u(x',t)\ 2 \(f)(y',z')\ 2 dr' 



e\r — r 



a 



Q[u]l 



[b,c] 



(3) 



The constant a is e 2 /C where C summarizes the result of the integration on the plane yz 
and an averaging of the potential over the width of the well. Dimensionally C is a length 
and defines the capacitance of the double barrier; 



QV 



\u(x, t)\ 2 dx 



(4) 



is the adimensional charge trapped in the well [6, c] at time t. Within the above approxima- 
tions we get the following one-dimensional nonlinear Schrodinger equation 



., d . , 
in—uix, t) 
dt v ' 



2m dx 2 



+ V(x) + aQ[u]l [bi , 



u(x } t) 



(5) 



which describes the dynamics of our quantum capacitor. An initial condition for u(x, 0) has 
to be given. We suppose u(x,0) to be a Gaussian wave packet moving toward the double 
barrier and initially centered far from the well [6, c] so that Q ~ at t = 0. 
Equation (5) has two conserved quantities, namely the number of particles 



N 



u(x, t)u(x, t) dx 



(6) 



and the energy 



E 



+ oo . 



U[X, t) 



2m dx 2 



+ V(x) 



u(x,t) dx + — aQ[w] 2 



(7) 



Some further comments are in order. The possibility of nonlinear effects in resonant 
tunneling was first envisaged in a paper by Ricco and Azbel [5] on the basis of a simple 
reaction-dissipation argument. The picture which emerges from [2] and the present study, 
as well as from a recent study of Malomed and Azbel [6] based on a different approach, is 
considerably more complex. 

The time scale of these effects is in principle within the reach of present experimental 
techniques. An experimental verification would be of special interest because as shown 
in [7] this type of concentrated nonlinearities can induce a genuine chaotic behavior in 
heterostructures which does not have an obvious classical counterpart. 

Equation (5) belongs to a class of equations which have other remarkable mathematical 
properties [8] beside those discussed in the present paper. 
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II. THE ONE-MODE APPROXIMATION 



In the following we will use the atomic unit system with h = 2m = 1. In that case 
lengths are given in Bohr radii, 1 a = 0.529 • 10~ 10 m, and energies in Rydbergs, 1 Ry= 13.6 
eV. The unit of time is 4.83 • 10" 17 s. 

We are interested in the initial value problem 

(A + D 2 X + V(x) + aQ[u] l [6 ,c]) u(x, t) = (8a) 



Q[ u ] = II u II [&,<=] = / 1^(^,01 dx (8b) 

' J b 

u(x,0) = fJe-i {x - Xo)2a ~ 2+tkoX (8c) 

where D t = -i§- v D x = and with V(x) = V (l[ a ,b] + l[c,d]) We consider a > 0, (3 > 0, 

a > 0, k > 0, a — x ^> <J and close to the real part of a resonance A(0) = Er(0) — iF(0) /2 
with < Er{0) < Vo and T(0) > 0, of the stationary Schrodinger operator D 2 X + V(x) on 1Z. 

In this section we shall describe a one-mode approximation to problem (8) which hope- 
fully can be rigorously justified in the limit <r, d — c — > oo when Vo and c — b are fixed 
and a, (3 are not too large, and provided that A(0) is a resonance so that Er(0) tends to 
some eigenvalue smaller than Vo of the potential Vo ^l]_oo,6] + l[c,oo[) an d T(0) — > when 
d — c — y oo. For simplicity we shall also assume that Er(0) is the only eigevalue of this 
potential at energy smaller than Vo- 

We observe that the problem (8) is formally a Schrodinger equation with a time de- 
pendent potential. The standard way to study this situation is to work with the in- 
stantaneous eigenstates and eigenvalues for the associated time dependent Hamiltonian 
H = D 2 X + V(x) + a<21 [6iC ] [9]. One says that the k-th. instantaneous eigenstate evolves 
adiabatically when transitions from it to any other ra-th eigenstate cannot occur. This 
happens when [9] 



[Eh — E, n 



< 1 (9) 



where Eh m is the matrix element of the operator d t H between the k-th and the ra-th eigen- 
states which have eigenvalues Eh and E m . In the problem (8) the instantaneous spectrum 
consists of a resonant state below the potential barrier Vo and a continuum band above Vo- 
We are interested to make a one-mode adiabatic approximation by keepeng only the in- 
stantaneous resonant state in the description of the wave function inside the double barrier. 
Since the distance between the resonant state and the rest of the spectrum is of the order 
of Vq, we get the adiabatic condition 



a 



Q 



< i • (io) 



We shall verify a posteriori that this condition is satisfied in the range of parameters we use. 
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We now derive the equation for the one-mode approximation. Let V(x) = Vo l[ a ,d] be 
the potential obtained from V(x) by filling the potential well [6, c\. We start by solving the 
initial value problem 

(D t + D 2 x + V(xj)ji(x,t) = (11a) 

jl(x } 0) = u(x } 0) (Hb) 
and look for a solution of Eq. (8) of the form u = jl + v. Then Eq. (8) becomes 

[D t + D 2 X + V(x) + aQ[Jji + v] l[ 6jC ]) v(x,t) = (V(x) - V(x) - aQ[Jj, + v] l[ 6jC ]) fi,(x,t) 

(12a) 

Q[J1 + v] = / \Ji(x, t) + v(x, t) | 2 dx (12b) 
u(a;,0) = 0. (12c) 



If < s < Vo, we let A(s) = Er(s) — iT(s)/2 be the resonance of D 2 x + V s (x) where 
V s (x) = V(x) + s l[b, c ]- Let e(s,x) be the corresponding resonance state 

(D 2 x + V s (x)-X(s))e(s,x) = . (13) 



We recall that e(s,x) is a function on 1Z such that e(s,x) oc exp yii^J X(s)xj , for x — > ±cx3. 
We also observe that, as a special case of the method of complex scaling [10], e(s,x) is of 
class L 2 on the contour 7 = e~ te ] — 00, a] \J [a,d] \J e 10 [0?, +oo[ if > is conveniently 
chosen. We normalize the resonance state to unity 



z(s,x) 2 dx = 1 . (14) 

In the spirit of the one-mode approximation, we assume v(x,t) = z(t)e(s, x) with 

s = olQ\Jjl + ?e(s)] . (15) 

The problem (12) then becomes 

[p t + D 2 X + V(x) + ctQ[fi + ze] 1[ 6)C ]) z(t)e(s, x) = (V - aQ[fi + ze\) l [6jC ] Jjl(x, t) (16) 

with initial condition z{0) = 0. If Eq. (15) defines a unique s = s[ji,z\, Eq. (16) reduces to 

(D t + X(s)) z(t)e(s, x) = (Vo - s) 1 M Ji{x, t) . (17) 

We multiply the above equation by e(s,x) and integrate with respect to x on the contour 
7. We also use the normalization condition for e implying / D t [e(s,x)] e(s,x) dx = and 
we obtain 
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(A + X(s))z(t) = (V - s) / p,(x,t)e(s,x) dx . (18) 

Jb 

Introducing = e tEot Jl(x,t) and = e tEot z(t), where Eq = k$, we can rewrite Eq. 

(15) and Eq. (18) as 

- -i (A(s) - E ) z(t) + i(V - s) I fj,(x,t)e(s,x) dx (19a) 



dt 

s = aQlfi + ze] (19b) 

*(0) = . (19c) 

The above three equations represent our one-mode approximation. Since s depends on zit) 
when a / 0, Eq. (19) is clearly a nonlinear ordinary differential equation with a driving 
term. 

A. The driving term 

In principle a complete integral representation of the driving term jl(x } t) or fi(x, t) could 
be given, but we prefer to use various approximations which seem justified in the limit <r, 
d — a — y oo. 

We start by looking for a solution of the stationary problem (^D^ + V(x) — E^j e(k, x) = 
of the form 

e tkx + r(k)e~ %kx x < a 

e(k, x) = —= < g + (k)e KX + g~(k)e~ KX a<x<d (20) 



'27T 



t(k)e lkx d < x 



where k = \[~E and k = \/Vo — k 2 . The normalization constant is chosen in such a way that 
/ e(k, x)e{k\ x) dx = S(k — k'). The functions r(k), t(k) and g ± (k) can be evaluated exactly 
by requiring e to be of class C l . In the limit of n(d — a) ^> 1 we get the simple semiclassical 
result 

ik + k 



r(k) = l -^le 2lka + O ( e -< d - a A (21) 
ik — k V / 



g -(k) = _^_ e («'*+«)« + O ( e -< d ~ a A (22) 



t(k), g + (k) = 0(e-« d -)) . (23) 
Now we have the solution of Eq. (11) of the form 

rco 

ji(x,t)= f{k)e~ lk 1 e(k,x) dk (24) 
Jo 
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where, in the limit k cr ^> 1, f(k) is the Fourier transform of /i(x,0) 

f(k) = p ae -h(^- k ?- 2 +^- k )^ . (25) 

The integral in Eq. (24) can be evaluated analytically by extending the lower limit to — oo 
and approximating the smooth functions of k with their value at k = k . The integration of 
the remaining quadratic exponential gives for a < x < d 

2^oe-"°(*-") (3 j [k a + i(a - xp)*- 1 ] 2 1, 22 ,., \ (0({ , 

U\x.t = , = exp < — k n cr + ik x > ZD 

PV ' ; ik -K Vl + 2ita- 2 F \ 2(l + 2zta- 2 ) 2 °J V ; 



where k = yVo — k$. This expression further simplifies for small times or large a. By 
expanding with respect to tcr~ 2 and using A; cr ^> 1 we get to the leading order 

H(x,t) = 1 7 e -2(-J 27 

zfc - K 

where t = (a — x )/(2k ) is the time a classical particle with kinetic energy E spends to 
cover the distance a — x and r = 2\/2k /(T is the full width at height e -1 / 2 of \f(k)\ 2 } i.e. 
the energy spread of the initial wave function u(x, 0). 

B. Solution of the one-mode approximation in the linear case 

We limit ourselves to a qualitative discussion of the limiting behaviors for r <C T(0) 
and r ^> r(0). Since a = 0, Eq. (19) becomes 

-- -i (MO) -E )z(t) + iV f C n(x,t)e(0,x) dx (28) 



dt 

with initial condition z{0) = 0. In the resonant case E = Er{0) we get 

+ ^^-z(t) = lV o F (t(x,t)e(0,x) dx . (29) 
dt 2 Jb 

The driving term in the right hand side has been evaluated in the previous section. It has 
a roughly costant argument and a modulus which behaves like exp | (r /2) 2 (t — t ) 2 ) at 
least for a large. In this case if r <C T(0) the Gaussian right hand side of Eq. (29) varies 
very little during the characteristic time period 2/T(0) of the left hand side and we get the 
approximate solution 

z{t) ^ f^j J b MM)e(M) dx (30) 

which has the same shape as the driving term, i.e. a Gaussian shape. On the other hand if 
r ^> r(0) we get the approximate solution 

z(t) ~ iVoe-^ T f C fi(x,t')e(0,x) dx dt' (31) 

JO Jb 
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characterized by a linear exponential decay slow with respect to the rise time. 
The above results can be generalized to the non resonant case. We get 

z(t) ~ — — — — — , , , / u,(x,t)e(0,x) dx (32) 

yj E R (0)-E(0)-ir(0)/2Jb PV J y ' ' 

when r <C T(O) and 

z(t) ~ iVoe-^ 11 ^-^-^)* f [ C ti(x,t')e(0,x) dx dt' (33) 

JO Jb 

when r > T(0). 

It is worth noting the conceptual simplicity of the one-mode approach in this case if 
compared with previous treatments [11]. 



C. Numerical simulations 

We now introduce the numerical simulations of the one-mode approximation obtained by 
a standard Runge-Kutta integration of Eq. (19). For comparison we give also the numerical 
simulations of the Eq. (8) obtained with the algorithm described in Appendix B. 

We consider an external potential V(x) with c — b = 15, d — c = b — a = 20 and Vo = 
0.3/13.6. In this case we have a resonace A(0) with Er(0) ~ 0.5 Vo and r — 5 • 10~ 3 /13.6. 
The initial wave function is chosen to be in resonance with = Er{0). We consider 
three different values for its width: a = 5775, a = 825 and a = 110 corresponding to 
T ~ 0.14 T(0), T ~ r(0) and r ~ 7.3 r(0), respectively. A constant (3 = (825^)" 1 / 2 
is used is every case. Finally we put a — x = 5<r so that no appreciable charge Q is 
present in the well at t = 0. Notice that the semiclassical condition and the condition 
k cr ^> 1, employed in evaluating the approximate expression (26) for the driving term, are 
well satisfied with the above parameters. 

Fig. 1, 2 and 3 show the behavior of the charge Q(t) for the three given values of a and 
for different values of the parameter a. Fig. 4, 5 and 6 show the trajectories of zit) in the 
complex plane for the same cases. We notice that the one-mode approximation reproduces 
exceptionally well the exact result. A small difference appears and grows up when a is 
decreased and a is increased in agreement with the adiabatic condition (10). This point will 
be discussed in detail in Section IV. 

The limiting behaviors for r <^ T(0) and r ^> T(0) in the linear case (a = 0) as 
expressed by Eq. (30) and Eq. (31) are apparent in Fig. 1 and Fig. 3. To understand qual- 
itatively the nonlinear behavior we introduce a simplified one-mode approximation which 
keeps the basic features of Eq. (19) and is analytically more manageable. 



III. THE SIMPLIFIED ONE-MODE APPROXIMATION 



If we assume that Eq. (19b) has a unique solution s = s[fi + ze], then the one-mode 
approximation (o.m.a.) defined by Eq. (19) can be written 

z = V(t,z(t)) (34) 
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where 

V(t, z) = —i (A(s[/i + ze]) — Eq) z + i(Vo — + ze]) I t)e(s[/i + ze], x) dx . (35) 



s=0 



x ^(0)cos 2 [y/EM0)(c-b)/2] 

y cos 2 [^(0)(c - b)/2] + E R (0)^V o - E R (0)(c - b)/2 



It turns out that the qualitative numerical aspects of the vector field V are already captured 
by a simplified one-mode approximation (s.o.m.a.). The s.o.m.a. is obtained from the o.m.a. 
i) replacing e(s, x) everywhere by e(0, x), ii) replacing Vo — s by Vo, in) replacing Q[ji -\- ze] 
by Qt- 26 ] an( i i v ) replacing A(s) — E by E' R (0)s — zT(0)/2, where the derivative 

is evaluated from Eq. (A4). The condition in Eq. (19b) now simplifies to 

s = aq\z\ 2 q = / |e(0, x) \ 2 dx (37) 

Jb 

and the vector field becomes 

V(t,z) = - + iE' R (0)aq\z\ 2 ) z + iV [ C fi(x,t)e(0,x) dx . (38) 



According to Eq. (27) during the interesting part of the evolution the driving term is well 
approximated by 

iV f C fx(x,t)e(0,x) dx = F e -K^) 2 (*-*°) 2 (39) 

Jb 

where 

f=« r e -M- Mo> „ & . (40) 

k — ik Jb 



A. solution of the simplified one-mode approximation 



We consider first two limiting behaviors in analogy to the discussion of the linear case 
in Section II B. 

When r <^ r(0) + E' R (0)aq\z\ 2 it is zit) ~ and zit) is close to the instantaneous fixed 
point z c (t) defined by V(t,z c (t)) = 0. If z c (t) is such a fixed point we have 



V 2 



fi(x } t)e(0, x) dx 



(41) 



and conversely if \z c \ solves Eq. (41) then we get the unique fixed point of modulus \z c 

iVo J 6 C fi(x, t)e(0, x) dx 
Z ° = r(0)/2 + iE' R (0)aq\z c \ 2 ' 



(42) 
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Since the left hand side of Eq. (41) is a strictly increasing function of \z c \ it is clear that V 
has a unique fixed point given by Eq. (41) and Eq. (42). In this situation we do not have 
oscillations in |z c (OI anc ^ therefore in Q(t) as shown in Fig. 1. 

Independently on the validity of the previous inequality, we distinguish two zones for the 
fixed point behavior, a) the linear zone : E R (0)aq\z c \ 2 <C T(0)/2. In this zone we have 

* c (f) ^ • (43) 

Using this expression the linear zone is defined by the condition 

8£^(0)ag|F|V(^) 2( *-* o)2 < T(0) 3 . (44) 

Due to the time dependence of the driving term we are sure to be in the linear zone at the 
beginning and at the end of the process, b) the nonlinear zone: E R (0)aq\z c \ 2 ^> T(0). In 
this zone from Eq. (41) we have 



Mt)| 

and from Eq. (42) 



Vq J 6 C (J,(x, t)e(0, x) dx l ^ 



E' R (0)aq 



(45) 



argz c (t) ~ arg / /i(x, t)e(0, x) dx . (46) 

Jb 

Therefore the nonlinear zone is defined by 

{E' R {0)aq) 1/3 |F| 2 / 3 e-K^)V*°) 2 » T(0) . (47) 



Notice that the argument of z c changes by tt/2 when we go from the linear to the nonlinear 
zone. When the driving term reaches its maximum, z c (t) becomes approximatively station- 
ary and then returns to the origin following a curve close to the outgoing trajectory. This 
kind of behavior is well illustrated in Fig. 4 and 5. 

When r ^> T(O) + E R (0)aq\z\ 2 we have the other limiting behavior characterized by the 
approximate equation 



z(t) ~ iVoe-^-^^fo W*')l 2 dt ' f f^^'je^j) dx dt 1 (48 



ft ,-c 
10 Jb 

and therefore 



\z(t)\ 2 = |F| 2 e- r (°)* QVK^)V-M 2 dt ^ (49) 

which is as in the linear case without oscillations in Q(t). The corresponding condition can 
be written more explicitely as 

r » T(0) + E' R (0)aq\F\ 2 e- r W QfVH^)V-*o) 2 dt ^ . (50) 
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The above condition is well satisfied in the simulation of Fig. 3 corresponding to a = 0.05 
while it is violated for a = 5. In this case oscillations of Q(t) appear. In the case a = 0.5 
the two sides of Eq. (50) are of the same order. From Fig. 6 we see that the trajectory of 
zit) is always far from the trajectory of z c (t). Still g|,2 c (t)| 2 gives a good estimate of the 
average, i.e. non oscillating, value of Q(t). 

For intermediate situations, that is r — T(0), the analytical treatment is necessarily 
more qualitative. In this connection it is convenient to introduce the linearized solution 
zi(t) around the fixed point z c (t) 




(51) 



where 




f I dz V c ' dz 



Qj{ Z c) d z { Z C 

The linearization matrix C has instantaneous eigenvalues 



(52) 



dV 

dz K ' 



which in the s.o.m.a. are given by 



\ 



dV , 

dz 



(53) 



£ ± = -E^&±iy/sE' R (0)aq\ 



(54) 



In the linear zone the real part dominates in this expression while in the nonlinear zone 
it is the imaginary part that dominates. In every case the fixed point is attractive since 
£± < 0. In the nonlinear zone we have an important rotation around the fixed point with 
angular speed 



S m £ ± = ^3 (E' R (0)aq) 1/3 \F 



1/3 | p,2/3 



(55) 



while the purely attractive behavior dominates in the linear zone. It is also possible to see 
that the rotation is oriented clock-wise. 

In Fig. 7 we give the result of the simulation of the s.o.m.a. for a = 825 and a = 5 
and we compare it with the linearization around the corresponding fixed point. We see a 
qualitative agreement with the o.m.a. result in Fig. 5. The linearization around the fixed 
point provides a good insight into the structure of the oscillating behavior. 



IV. ADIABATIC CONDITION 

In this section we give an explicit expression of the adiabatic condition (10) in terms 
of the relevant parameters. We use the s.o.m.a. with Q(t) = q\z(t)\ 2 and we situate 
our discussion in the most interesting region where \t — t \T /2 < 1. We approximate 
d\z(t) \ 2 / dt ~ d\z c (t) \ 2 /dt. This is not entirely correct because, due to the oscillations, the 
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derivative of |^(t)| 2 may reach higher values. However, and this is confirmed by numerical 
simulations, the order of magnitude of the maximum value of d\z(t) \ 2 / dt is given by the slope 
of the first oscillation which is close to the maximum value of d\z c (t) \ 2 /dt. In the interesting 
nonlinear zone the value of |z c (t)| 2 is given by Eq. (45) and its time derivative is evaluated 
by using Eq. (39) for the driving term 



d\z c {tW 



dt 



< r 



I ^ 1 2/3 



a 



2/3 



(56) 



where we used E' R (0)q ~ 1. Noticing that e(0,x) ~ (c — b)~ l l 2 in [6,c], from Eq. (40) we get 



F 



K Q \fc 



(57) 



The adiabatic condition (10) becomes 

1/3 r (2f3k ^V e-^ 



(b-a) \ 2 / 3 



a 



< 1 



(58) 



We make this inequality physically more transparent by inserting the width of the resonance. 
From Eq. (A9) we have 



r(o) 



K 2 e -2K (b-a) 



k (c - b) 

and, since in our case k ~ /c , the adiabatic condition (10) is also 

1/3 



(59) 




V 




< 1 



(60) 



The order of magnitude of the left hand side is in good agreement with the maximal value 



of a 



Q V as obtained from numerical simulations and the inequality is well satisfied. In 



fact the left hand side ranges from 10 in the case a = 5775 to 3 • 10 in the case a = 110 
for a = 5. 
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APPENDIX A: EIGENVALUE AND EIGENFUNCTION AT RESONANCE 

We want to calculate the eigenvalue A(s) = Er(s) — iT(s)/2 and the eigenfunction e(s, x) 
of the stationary problem 



12 



(Dl + V a (x)-X(s))e(s,x) = (Al) 

with V s (x) = Vo l[ a ,b] + s l[6,c] + Vo l[c,d] and < s < Vq. The eigenfunction e(s, x) has to be 
known only in the interval [6, c\. For convenience we choose a = —d and b = —c. 

Working with the assumption that d — c is large, we approximate Er(s) by the first 
eigenvalue E for the confining potential V SjO0 (x) = Vo l]-oo,&] + s l[6,c] + Vo l[c,oo[ an d we 
approximate the resonant state e(s, x) in [6, c] by the corresponding eigenfunction. Observing 
that e(s,x) is an even function of x, up to a normalization constant we have 

, , f cos(a/ E — s x) < x < c . . n . 

The C 1 requirement at i = c gives 

7 = cos - s c) (A3) 

and 



V )^lf = tan (VE^~s c) . (A4) 
\JE — s v 1 

This equation has a unique solution s < E < Vo which can be found numerically. Once the 
solution E is known, we normalize the corresponding eigenfunction and we get 

cos(\/ E — s x) , , , 

e( s ,x) = — j = y ' = (A5) 

, sin(2V_E-s c) . l+cos(2V-B-s c) 



2VS :: ^ 2VVo--B 

in the interesting interval < x < c. 

Finally we derive a semiclassical formula for Sm A(s). In this case we improve the above 
approximation for e(s,x) valid only for \x\ < c by considering the reflected waves in the 
potential V s (x). In the semiclassical limit up to a normalization constant we can write 

cos(v / E — s x) < x < c 

x) = \ 7 e-^^ (*" c ) + Se-^^v (d-c)+Vv^E ( x -d) c<x<d . ( A6 ) 

^iy/E (x-d) d < X 

The costant 7 is given by Eq. (A3) while 8 and rj can be obtained by C 1 requirement at 
x = d. In particular we have 

V = cos ( y^ c) e"^ . (A7) 

V Vo - E - i\/ E 

By using Eq. (Al) and its complex conjugate we can write 
[A(s) — A(s)j y _ e(s, x) e(s, x) dx = 

/d rd r _ — " 

„e(s,x) e(s,x) dx — j _d^.e(s,x) e(s,x) o?x = e(s,d) d x e(s,d) (A8) 

-ci J — ci L J 
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where d — d > is small. Inserting the expression (A6) for e(s, x) at x = d we get 



see 



-2y/%=E (d-c) 



t/ I i sin(2\/E — s c) . l+cos(2\/E — s c) 
V ° l C+ 2^ + 



2-yVo^E 



(A9) 



where _E = Er(s) has been evaluated above. 



APPENDIX B: NUMERICAL ALGORITHM 

Here we discuss briefly the algorithm used for the numerical simulations of Eq. (8). The 
partial differential equation is transformed in a set of algebraic equations by discretizing the 
space-time x — t in an appropriate two dimensional lattice. Let us suppose that we want to 
follow the temporal evolution of the initial state u(x, 0) during the time interval [0, t max \. Let 
[xmim %max] be the space interval where the wavefunction u(x, t) can be considered localized 
for < t < t max . The space-time lattice is defined by putting 



x x mm + j Ax j = 0, 1, 2, . . . , J + 1 



fBl] 



t^nAt n = 0,1,2,..., N 



(B2) 



where the total number of points in the lattice is related to the lattice constants Ax and At 
by 



(B3) 



The wave function reduces on the lattice to a matrix 



(B4) 



u[x, t) — y u n 3 



(B5) 



where u° are kown. The Hamiltonian operator H = D 2 X + 14 7 , where W = V + ctQl[b, c ] 
includes both the external and the nonlinear potentials, becomes 



Hu(x,t) {Hu) n 3 = -Ax- 2 {u a J+l - 2u] + u]^) + Wfu, 



(B6) 



with Wf = W(x h t n ). 

The discretization of the operator D t is obtained with the Cayley approximation for the 
exponential time evolution operator 



-iAtH ^ 



1 + -AtH 
2 



1 - -AtH 
2 



(B7) 



which preserves unitarity. The finite difference scheme for time evolution gives the implicit 
equation 
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1 + -AtH 

2 



u(t + At) 



1 - -AtH 

2 



«(*) 



fB8) 



which corresponds in the lattice to the following linear system 



uu 



n+l 



uu 



n+l 
J + l 



(B9) 



where 



Xi 



uu- 



;bio) 



anc 



At 



A; 



;2 + g AfW? 



At 



2 A; 



fB12) 



By using the condition u r 



u 



localized in \x 



mm •} ^max] 



J+l 



which expresses the fact that the wavefunction is 



we can write the following tridiagonal system 





u 


. 


.. 








u 


tn 
$2 


u . 


.. 











U 


tn 

?3 • 


.. 














. 


tn 
■ ■ Cj-2 


u 











. 


u 


tn 

Cj-i 


u 


V o 





. 


.. 


u 


tn 



/ u l +1 \ 

«2 + 1 
, n+l 



n+l 
U J-2 



( x n i \ 

Xs 

n 

Xj-2 

xU 

\ X n j / 



;bi3) 



Starting from u° the solution of the above system at each time step gives the unknown it" +1 
in terms of the known it". 

Due to its tridiagonal nature the linear system in Eq. (B13) is efficently solved by stan- 
dard triangularization methods requiring a number of operation O(J) with a total compu- 
tation time O(JN). The number of lattice points N and J is chosen in such a way that 
it™ is a good approximation to the continuous function u(x,t). With the parameters given 
in Section II, the difference between it" and u(x,t) is estimated to be less than 0.1 % for 
Ax < 1 and At < 4. 
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FIGURES 



FIG. 1. Comparison between the charge Q(t) obtained by numerical solution of the partial 
differential equation (8) (dot-dashed line) and of the one-mode approximation (19) (full line) as a 
function of the nonlinearity strength a. Notice that the two curves almost coincide. The values of 
the parameters are given in the text. Here a = 5775 which corresponds to r — 0.14 T(0). 

FIG. 2. As in Fig. 1 but for a = 825 which corresponds to r — T(0). 

FIG. 3. As in Fig. 1 but for a = 110 which corresponds to r — 7.3 T(0). 

FIG. 4. Trajectory of the numerical solution z(t) of the one-mode approximation (19) for the 
same cases of Fig. 1 (full line). The dashed line is the trajectory of the corresponding fixed point 
z c (t). 

FIG. 5. Trajectory of the numerical solution z(t) of the one-mode approximation (19) for the 
same cases of Fig. 2 (full line). The dashed line is the trajectory of the corresponding fixed point 
z c (t). 

FIG. 6. Trajectory of the numerical solution z(t) of the one-mode approximation (19) for the 
same cases of Fig. 3 (full line). The dashed line is the trajectory of the corresponding fixed point 
z c (t). 

FIG. 7. Trajectory of the numerical solution z(t) of the simplified one-mode approximation for 
the case a = 825, a = 5 (full line). The dashed line is the trajectory of the corresponding fixed 
point z c (t) while the dotdashed line is the linearized solution Zi(t). 
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